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ABSTRACT 



We analyse with simple real-space statistics the Virgo consortium's 
cosmological N-body simulations. Significant clustering rapidly develops well 
below the initial mean interparticle separation (Aj), where the gravitational force 
on a particle is dominated by that with its nearest neighbours. A power-law 
behaviour in the two point correlation function emerges, which in the subsequent 
evolution is continuously amplified and shifted to larger scales, in a roughly 
self-similar manner. We conclude that the density fluctuations at the smallest 
scales due to the particle-like nature of the distribution being evolved are thus 
essential in the development of these correlations, and not solely, as usually 
supposed, the very small continuous (fluid-like) fluctuations at scales larger than 
(A,,). 
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One of the central projects of contemporary cosmology is to explain within a single 
framework the small fluctuations observed in the temperature of the cosmic microwave 
background radiation (CMBR) and the large fluctuations observed in the distribution of 
galaxies and other discrete objects. The current paradigm ( |Peebles, 1993 , Fadmanabhan 



|1993|) for this framework is the cold dark matter (CDM) cosmology, in which it is 
(essentially) the growth under self-gravity of very tiny initial perturbations in a collisionless 
fluid which give rise to both kind of phenomena. The CMBR fluctuations are explained - 
with apparently great observational success - by the very tractable evolution of the fluid 
equations linearized in the perturbations. The extension to the formation of objects and 
their distributions is much more difficult, involving the highly non- linear regime of the 
equations. In practice the main theoretical instrument in this regime is numerical N-body 
simulation (NBS) ( |Hockney and Eastwood, 1981| ). From some given initial conditions (IC 
- generally a Gaussian random field with small-amplitude correlated fluctuations) specified 
by the linear evolution (and thus in principle constrained by observations of the CMBR) 
these NBS are supposed to model the evolution of the continuous density field under 
gravity into the non-linear era. Ultimately catalogues of different objects are extracted for 
comparison with observations, and the idea is that significant information about the IC 
in the continuous density field can be extracted by such a comparison. In this letter we 
re-examine the dynamics in the non-linear regime of some of the most elaborate cosmological 
NBS performed to date, those of the Virgo consortium ( | Jenkins et al., 19"98| ). We conclude 



that the correlations which emerge at these scales are not correctly understood solely in 
terms of the amplification of the initial small fluctuations at large scales, but rather involve 
essentially the large fluctuations, associated with the discreteness of the simulated system, 
at the very smallest scales. We discuss the important consequences of this observation. 



There are several widely used types of cosmological NBS ( |Hockney and Eastwood, 



1981| ): 'particle-mesh' (PM), 'particle-particle, particle-mesh' (P 3 M), 'adaptive mesh 
refinement' and tree codes. The Virgo consortium's NBS are of the adaptive P 3 M type, and 
our discussion here applies directly only to this type of NBS. The difference between the 
various methods has to do with how the gravitational force acting on a particle is modelled. 
The PM type computes the force from a potential calculated on a mesh, typically of the 
same size as the initial mean interparticle separation, while the P 3 M and other methods 
add to this the direct contribution to the force from particles on scales smaller than this 
distance. With either algorithm the relation between the simulated system and the physical 
system which one is attempting to model - CDM particles evolving under their self-gravity - 
is not at all evident. In NBS the mass of the 'particles' is (very) macroscopic, typically of the 
order of a galaxy mass, and not the microscopic mass of a CDM particle. PM NBS attempt 
to reproduce a treatment of the CDM as a continuous fluid, with the 'particle' putatively 
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representing an element of the fluid. The other methods of NBS effectively describe the 
dynamics of point particles interacting through gravity in an expanding Universe. In 
particular for P 3 M the 'particle-mesh' part in the sum of the gravitational force can be 
understood simply as a sufficiently accurate and efficient numerical way of calculating the 
force due to points outside the mesh cell, inside of which the force is calculated exactly 
with a smoothed 1/r potential. This latter smoothing can be understood as an effective 
'smearing' of the particles like in PM, but now at a characteristic scale e which, as we will 
see further below, is always much smaller than the mesh size. Physical results are assumed 
to be those independent of this smoothing scale, and in fact it can be understood as a 
purely numerical artefact introduced to avoid the divergence of the numerical time step in 
the integration of the equations of motions. 



In the cosmological context P M algorithms have been widely used ( JJenkins et al., 



1998| , Pavis et al., 1985j |Efsthathiou et al., 19871) because they allow a much greater range 



of spatial resolution while remaining numerically feasible. Despite its evident importance 
the central issue of how this or indeed the PM type of NBS describe CDM has received 
surprisingly little attention in the literature. Those few studies which have been undertaken 
of some of these issues ( [Melott, 1990| , [Khulman, Melott fc Shandarin, 1986| , (Sprinter et al., 



|1998| ) clearly support the idea that discreteness can be very important indeed (although 
possibly less so in PM simulations with a mesh size well above the interparticle separation). 
Our aim here is not to try to resolve this general question. Rather we are interested in 
understanding what the dynamics is which gives rise to the correlated structure in the 
non-linear regime which is observed to emerge in P 3 M NBS, which as we have remarked 
are simply simulations of self-gravitating particles subject to some particular initial 
perturbations. 

Usually the main instrument for characterising the evolution of structure in NBS is 
the power spectrum P(k) = \S(k)\ 2 where 6(k) is the Fourier transform (FT) of the density 
contrast field. In NBS it is calculated by assigning a local mass density on a grid, and 
then performing a fast FT. Here our approach is to look directly at real space quantities, 
calculated directly from the unsmoothed particle distributions. Given what we have said 
about the nature of the P 3 M NBS, this is a reasonable, and, as we will see, instructive 
approach. It is the one used in a recent study flBottaccio et al., 2002|) through NBS of 



gravitational dynamics in a purely Newtonian (non-expanding) system of particles with 
Poissonian IC. In this case the method has been used to establish the essential role played 
by particle fluctuations in the formation of structures. 



We have analyzed three cosmological NBS performed by the Virgo consortium flJenkins 



et al., 1998|) . They are all CDM models with N = 256 3 particles and gravitational force 
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smoothing length e = 0.036 Mpc/h 0. They all have the same value of the cosmological 
parameters, with in particular total mass density £1 = 1. The first two NBS, which we 
refer to as S-I and S-II, are of the "standard" CDM type, differing only in their resolution 
(number density): For S-I the volume V of the simulation is a cube of side L = 239.5 Mpc/ h, 
while for S-II L = 84.5 Mpc/ h. The third NBS, S-III, is the same as S-I, except that its 
IC are (slightly) different as they are those prescribed by a rCDM model. In S-I and 
S-III each particle then represents a mass of m p = 2.27 x 10 11 solar masses while in S-II 
trip = 1.0 x 10 10 . In S-I and S-III we find that the initial mean interparticle separation 
(A;) « 1 Mpc/h, while in S-II (A;) « 0.33 Mpc/h. Hence (A,) > e in all three NBS as 
noted above. The NBS are all run from a red-shift of z = 50 until today (z = 0), i.e. for a 
time which is essentially the age of the Universe I. This means that the time of evolution is 
essentially one dynamical time r dyn wl/ \fGp~, where the latter is simply the characteristic 
time scale associated with a mass density p = Nm p /V. A particle typically moves by a 
distance of order the lattice spacing in this time. 



Let us first consider S-I. The first statistic we consider is the conditional density ( Bylos 



[Labini, Montuori and Pietronero, 1998 , Gabriclli, Joyce fc Sylos Labini, 2002 ) given by 



r(r) = (n(r)n(O)) / (n) where n(r) is the microscopic number density. It is simply the mean 
density of points at a distance r from an occupied point. It is plotted in Fig.|l| for a sequence 
of time slices, beginning from the initial distribution at z = 50 until today. Also shown is 
the pre-initial 'glass' configuration. This pre-initial distribution (often taken as a regular 
lattice in NBS) is obtained by evolving the simulation with the sign of gravity reversed, 
which makes the particles settle down to a "uniform" configuration in which the force on 
each particle is very close to zero flJenkins et al., 1998j ). This distribution is very isotropic 
but it is still characterized by long-range order similar to that of a regular lattice (|Gabrielli,| 
Joyce fc Sylos Labini, 2002|) . It has correlated fluctuations at all scales, and in particular 



large fluctuations (of order one) at scales of order (Aj). The IC of the NBS at the initial 
redshift z = 50 are obtained by applying small (compared to (Aj)) displacements to the 
points in this pre-initial distribution, the displacement field being prescribed by the power 
spectrum of the (continuous) theoretical model i. The behaviour of T(r) for the glassy 



1 With this smoothing the force is 53.6% of the true 1/r 2 force at e and more than 99% at 2e (Jenkins et 



al., 1998[ ). 

2 In these matter dominated cosmologies the time is given by t = t j (1 + z) 3 / 2 where t a is the age of the 
Universe today. 

3 A real space analysis of these IC complementary to that here is given in ( Baertschiger & Sylos Labini 



2002). Very significant differences between these and those expected from the theoretical input power 
spectrum are found even at large scales. 
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configuration and these IC at z = 50 are very similar. This is because the perturbations are 
very small in amplitude compared to the initial interparticle separation (Aj). The highly 
ordered lattice- like nature of the initial distribution is manifest: there is an excluded volume 
around each point so that T(r) is negligible until very close to (Aj), where it shows a small 
peak, with some oscillation about the mean density evident corresponding to the long-range 
order (|Gabrielli, Joyce fc Sylos Labini, 2002|) . When the evolution under self-gravity starts 
the excluded volume feature is rapidly diluted, having almost completely disappeared by 
z = 3 (t = 0.125t o ). This corresponds to the development of power-law clustering at very 
small scales where it was completely absent in the IC. 

Let us consider what is driving this dynamics. In Fig.0 we show, for the pre-initial and 
initial perturbed configuration, the mean value of the magnitude of the gravitational force 
on a given point due to the mass inside a sphere of radius R centred on it. In the pre-initial 
distribution the force is indeed zero within the errors, as expected. In the perturbed 
configurations one sees that the "long-range" contribution dominates only marginally (by 
a factor of about three) over that coming from the nearest neighbour (NN) ( |Baertschigcrj 
fe Sylos Labini, 2002|) . On the figure is shown also the force on a point due to a single 
other one at a distance (Aj): it is, as expected, of order the force due to the nearest points. 
When points move towards their NN the first contribution grows as one over the distance 
squared while the latter changes, supposing that the validity of linear growth of fluctuations 
at larger scales, in proportion to the scale factor, i.e. as 1/(1 + z). Thus at z = 5 we would 
expect the latter to grow by a factor of 10 and thus the NN force to dominate below about 
(Aj)/5. That this is indeed the case is very clearly confirmed by considering the evolution 
of the NN distribution during the simulation. In the insert panel of Fig.|2] the probability 
that the NN of a particle lies at a distance r is shown, again in different time slices. The 
initial distribution (not shown) is extremely peaked about (Aj), then broadens and develops 
a second peak at small scales. This second peak at small scales is a clear "fingerprint" for 
the dominance by NN interactions: it is approximately at the scale 2e where the force is 
strongly cut off by the smoothing. The build-up of neighbours around this scale is due to 
the fact that the force which is driving their motion feels the lower cut-off in it at this scale. 

Between z = 5 and z = 3 we see - in this range of scales completely dominated by 
NN interactions - the appearance of an approximately power-law behaviour in T(r). At 
z = 3 the amplitude of T(r) over most of the range r < (Aj) is greater than the mean 
density so that this power law is also now seen in the normalised correlation function 
£(r) = r(r)/ (n) — 1 shown in the inset of Fig. 1. At the next time slice, at z — 1, the form 
of T(r) up to approximately 0.3 (Aj) is almost precisely the same, being simply amplified 
by an order of magnitude. At the same time the power law extends slightly further before 
flattening from V pa 3{n) to a smooth interpolation towards its asymptotic value. The 



- 6- 



evolution in the remaining two time slices is well described over the range T(r) > 3(n) as 
a simple amplification, and translation towards larger scales, of T(r). At the final time the 
power law (with exponent 7 ~ —1.7) extends to ~ 2(Aj). 

Let us now also consider the other NBS: S-II and S-III. In both cases we find the same 
qualitative evolution (see Fig.^J) as observed in the S-I simulation in the range of scales we 
are discussing. The principal difference is in the dependence of the amplitude and extent of 
the power-law type behaviour in T(r) as a function of time. In S-II the clustering develops 
more rapidly initially at small scales, while S-III lags behind relative to S-I. The reason for 
this difference is simply traceable to the relative amplitude of the perturbations in each 
case to the pre-initial distribution: S-II has the same perturbation field superimposed on a 
distribution with a smaller initial interparticle separation 0.33(Aj), while the perturbation 
field at the scale of the lattice in S-III is smaller in amplitude than in S-I. Thus it takes 
shorter (or longer) for the NN contribution to the force to dominate and the evolution of 
clustering due to these forces starts faster or slower. 

Apart from this difference the evolution is however strikingly similar. Note, for 
example, the almost perfect superposition of T(r) for S-III at z = and S-I at z — 0.3. The 
evolution of clustering in the non-linear regime thus appears to be almost independent of 
the IC, except trivially as represented by the difference in amplitude. Some dependence on 
the IC is evident at larger scales, but it is a weak residual dependence imposing a small 
deviation from the power-law which is formed at and then translated from the smallest 
resolved scales. In Bottaccio et al. (2002) a very similar behaviour is observed in the 
evolution of clustering of particles by Newtonian forces, without expansion and with simple 
Poissonian IC. The authors give a physical interpretation of this clustering in terms of 
the exportation of the initial "granularity" in the distribution to larger scales through 
clustering. The self-similarity in time of T(r) is explained as due to a coarse-graining 
performed by the dynamics: one supposes a NN dynamics between particle-like discrete 
masses, with the mass and physical scale changing as a function of time as the clustering 
evolves (particles forming clusters, clusters forming clusters of clusters etc.). This would 
appear to provide a good explanation for the behaviour observed here as well, but further 
theoretical work is clearly needed to establish this and find a more quantitative description. 

The primary conclusion we draw from our real space analysis of the Virgo NBS is 
therefore that the fluctuations at the smallest scales in these NBS - i.e. those associated 
with the discreteness of the particles - play a central role in the dynamics of clustering in 
the non-linear regime. In particular the power-law type correlations appear to be built up 
from the initial clustering at the smallest scales. The nature of the clustering (in particular 
the exponent of the power-law) seems to be independent (or at most very weakly dependent 
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on) the IC and its physical origin should therefore be explained by the dynamics of discrete 
self-gravitating systems. The fluid-like statistical description ( |Gabrielli, Joyce fc Sylos 



Labini, 2002[) and equation of motions, which is the framework used to describe a CDM 



universe, do not consider the non-analytical "particle" term of noise which is represented 
by NN interactions. The latter as we have seen are strongly present in the NBS we have 
analysed and appear to play a crucial role in the formation of the correlated structures 
observed to emerge. It is instructive to note that, in the paradigmatic example of stochastic 
(homogeneous and isotropic) point processes, the Poisson distribution, the gravitational 
force on an average point is due very predominantly to by its NN ( |Chandrasekhar, 1943|) . 



Large-scale small-amplitude density fluctuations do not give rise to a significant contribution 
to the force acting a point, because of symmetry: i.e. large-scale isotropy ( |Chandrasekhar 

rem 



This conclusion is in disagreement with the usual interpretation given of the origin 
of these correlations: they are supposed to represent solely the evolution in the non-linear 
regime described by a set of fluid equations, with the small fluctuations of the theoretical 
input model as IC. The dynamics we have described is essentially dependent on the 
gravitational forces at the smallest resolved scale in the NBS, and the small smooth 
component of the force added to this by the perturbations at larger scales appears to 
be irrelevant in the non-linear regime. What dominates the dynamics are the large 
fluctuations at small scales intrinsic to the discrete pre-initial distribution, and the imposed 
fluctuations play only a minor role in the non-linear regime where the fluctuations are 
large. In cosmology the link between the IC and the "predictions" for structure formation 
in the non-linear regime is very important, as it is through it that one tries to constrain 
models using observations of galaxy distributions. Our conclusions completely change the 
perspective on this problem. Power-law correlations of this type are the most striking and 
well established feature of such distributions ( [Peebles, 1993| , [Sylos Labini, Montuori and 



[Pietronero, 1998| ). The theoretical problem of their origin therefore must deal with the 
apparently crucial role in their formation of an intrinsically highly fluctuating (and thus 
non-analytic) density field. In particular the origin of the exponent in the correlations 
and the dependence of the extent of such correlations on the discretization (physical or 
numerical) needs to be understood. In particular we note that, despite the centrality of 
discreteness, the clustering found in the Virgo NBS (and those of |Bottaccio et al., 2002|) is 
almost independent of the smoothing scale e which cuts off the force at small distances. We 
will be studying this and other issues in future work. 

We thank M. Bottaccio, H. De Vega, R. Durrer, A. Gabrielli, A. Melott, M. Montuori, 
L. Pietronero and N. Sanchez for useful discussions. F.S.L. acknowledge the support of the 
Swiss National Science Foundation. 
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Fig. 1. — Two-point conditional density for the NBS S-I. The initial mean interparticle 
separation (Aj) and the softening length e are shown, as well together with the best power- 
law fit at the end of the simulation. In the insert panel the evolution of the reduced 
correlation function £(r) is shown by vertical dashed and dashed-dotted lines respectively. 
The arrows indicate qualitatively how clustering proceeds: first it develops at small scales as 
some particles move towards their nearest neighbor, and then it is subsequentely exported 
to larger scales. 

Fig. 2. — The behaviour of the average total force on a point (and its variance) due to the 
mass in a sphere of given radius r about it, for the pre- initial (glass) and initial particle 
distribution. The contribution from a single particle at distance (Aj) is also shown. The 
units are given by the choice Grn^ — 1. In the insert panel it is shown the evolution of the 
nearest neighbor distribution in S-I. The vertical dashed marks (Aj). 

Fig. 3. — Comparison of the time evolution of T(r) in S-I, S-II and S-III, normalized to 
the mass density. For S-II and S-III there are shown (in closkwise order) the time slices for 
z = 10, 1, while for S-I it is also shown the slice at z — 0.3. (A) LR is the average distance 
between nearest neighbors in the low-resolution simulation (S-I and S-III), while {A) HR is 
that scale in the high- resolution simulation (S-II). 
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